In this script we are going to map T-cell subtypes onto the Visium slides.
library(Seurat)
library(ggpubr)
library(cowplot)
library(dplyr)
library(ggplot2)
library(stringr)
library(readr)
library(SPOTlight)
library(SPATA2)
Loading necessary paths and parameters
source(here::here("misc/paths.R"))
source(here::here("utils/bin.R"))
"{cd4}/{plt_dir}" %>%
glue::glue() %>%
here::here() %>%
dir.create(path = ,
showWarnings = FALSE,
recursive = TRUE)
"{cd4}/{robj_dir}" %>%
glue::glue() %>%
here::here() %>%
dir.create(path = ,
showWarnings = FALSE,
recursive = TRUE)
set.seed(123)
Extract sample id and get Donor ID
# sample_id <- params$sample_id
sample_id <- "esvq52_nluss5"
donor_id <- id_sp_df[id_sp_df$gem_id == sample_id, ] %>% dplyr::pull(donor_id)
Set gene dictionary
"{cd4}/gene_dict.R" %>%
glue::glue() %>%
here::here() %>%
source(file = .)
We have 8 different datasets that we are going to analyze separately. The spatial data comes from the script 03-clustering/03-clustering_integration.Rmd while the sc data can be found in Ramon’s scRNAseq analysis: /scratch/devel/rmassoni/tonsil_atlas_private/2-DOWNSTREAM_PROCESSING/results/R_objects/processed_seurat_objects/processed_seurat_objects/tonsil_integrated_with_harmony_scrublet_annotated.rds.
sp_obj <- "{clust}/{robj_dir}/integrated_spatial_annotated.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)
# Load SPOTlight data
spotlight_ls <- "{cd4}/{robj_dir}/spotlight_ls_cd4.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)
Load MAGIC data from the script MAGIC_denoising.Rmd
magic_df <- "{cd4}/{robj_dir}/MAGIC-mtrx.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)
# create a new assay to store ADT information
magic_assay <- CreateAssayObject(counts = as.matrix(magic_df))
# add this assay to the previously created Seurat object
sp_obj[["MAGIC_Spatial"]] <- magic_assay
Seurat::DefaultAssay(sp_obj) <- "MAGIC_Spatial"
Add colors to cell types
col_vec <- c(
"Naive" = "gray88",
"CM Pre-non-Tfh" = "gray71",
"CM PreTfh" = "gray59",
"Tfh T:B border" = "#67a9cf",
"Tfh-LZ-GC" = "#3690c0",
"GC-Tfh-SAP" = "#02818a",
"GC-Tfh-0X40"= "#016c59",
"Tfh-Mem" = "#014636",
"T-Trans-Mem" = "#fd8d3c",
"T-Eff-Mem" = "#e31a1c",
"T-helper" = "#800026",
"Eff-Tregs" = "#df65b0",
"non-GC-Tf-regs" = "#e7298a",
"GC-Tf-regs" = "#ce1256"
)
(nm_df <- data.frame(col_vec))
## col_vec
## Naive gray88
## CM Pre-non-Tfh gray71
## CM PreTfh gray59
## Tfh T:B border #67a9cf
## Tfh-LZ-GC #3690c0
## GC-Tfh-SAP #02818a
## GC-Tfh-0X40 #016c59
## Tfh-Mem #014636
## T-Trans-Mem #fd8d3c
## T-Eff-Mem #e31a1c
## T-helper #800026
## Eff-Tregs #df65b0
## non-GC-Tf-regs #e7298a
## GC-Tf-regs #ce1256
nm_df$plt_nm <- rownames(nm_df)
nm_df$mod_nm <- stringr::str_replace_all(
string = rownames(nm_df),
pattern = "[[:punct:]]|[[:space:]]",
replacement = ".")
decon_mtrx <- spotlight_ls[[2]]
decon_mtrx <- decon_mtrx[, colnames(decon_mtrx) != "res_ss"]
# Set as 0 cell types predicted to be under 1 % of the spot
# decon_mtrx[decon_mtrx < 0.03] <- 0
Change column names
new_cn <- data.frame(mod_nm = colnames(decon_mtrx)) %>%
dplyr::left_join(nm_df, by = "mod_nm") %>%
dplyr::mutate(plt_nm = dplyr::if_else(is.na(plt_nm), mod_nm, plt_nm)) %>%
dplyr::distinct() %>%
dplyr::pull(plt_nm)
colnames(decon_mtrx) <- new_cn
We are going to add the deconvolution to the Seurat object.
sp_obj@meta.data <- cbind(sp_obj@meta.data, decon_mtrx)
Subset sample of interest
sp_sub <- sp_obj[, sp_obj@meta.data$gem_id == sample_id]
sp_sub@images <- sp_sub@images[Seurat::Images(sp_sub) == sample_id]
Check Topic profiles
nmf_mod_ls <- spotlight_ls[[1]]
nmf_mod <- nmf_mod_ls[[1]]
h <- NMF::coef(nmf_mod)
rownames(h) <- paste("Topic", 1:nrow(h), sep = "_")
topic_profile_plts <- SPOTlight::dot_plot_profiles_fun(
h = h,
train_cell_clust = nmf_mod_ls[[2]])
topic_profile_plts[[2]] +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90),
axis.text = ggplot2::element_text(size = 12))
Look at all cells profiles
topic_profile_plts[[1]] +
ggplot2::theme(
axis.text.y = ggplot2::element_blank(),
axis.text.x = ggplot2::element_blank(),
axis.title = ggplot2::element_blank())
Look at cells topic profile
basis_spotlight <- data.frame(NMF::basis(spotlight_ls[[1]][[1]]))
train_labs <- spotlight_ls[[1]][[2]]
colnames(basis_spotlight) <- unique(stringr::str_wrap(train_labs, width = 30))
basis_spotlight[basis_spotlight < 0.0000001] <- 0
DT::datatable(basis_spotlight, filter = "top")
Look at the location of each cell type in each slice separately
# Iterate over cell types
ct <- colnames(decon_mtrx)
# Iterate over images
lapply(names(sp_obj@images), function(nm) {
print(nm)
nm_donor <- id_sp_df %>% dplyr::filter(gem_id == nm) %>% dplyr::pull(donor_id)
# Iterate over cell types
ct_plt_ls <- lapply(ct, function(i) {
tmp_plt <- Seurat::SpatialFeaturePlot(
object = sp_obj,
features = i,
alpha = c(0, 1),
images = nm) +
ggplot2::scale_fill_gradientn(
colors = heat.colors(10, rev = TRUE)) +
ggplot2::scale_alpha(range = c(0, 1)) +
ggplot2::labs(title = stringr::str_wrap(string = i,
width = 25),
fill = "") +
ggplot2::theme(
plot.title = ggplot2::element_text(
hjust = 0.5,
size = 20,
face = "bold"))
if (sum(sp_sub@meta.data[, i]) == 0) {
tmp_plt <- suppressMessages(tmp_plt + ggplot2::scale_alpha(range = c(0,0)))
}
return(tmp_plt)
})
(plt_arr <- cowplot::plot_grid(
plotlist = ct_plt_ls,
axis = "trbl",
align = "hv",
nrow = 5,
ncol = 5))
"{cd4}/{plt_dir}/cell_type_location_cd4_{nm_donor}.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = plt_arr,
base_height = 25,
base_width = 25)
})
## [1] "tarwe1_xott6q"
## [1] "c28w2r_7jne4i"
## [1] "esvq52_nluss5"
## [1] "p7hv1g_tjgmyj"
## [1] "gcyl7c_cec61b"
## [1] "zrt7gl_lhyyar"
## [1] "qvwc8t_2vsr67"
## [1] "exvyh1_66caqq"
## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/CD4-Analysis/2020-09-22/plots_2020-09-22/cell_type_location_cd4_BCLL-2-T.pdf"
##
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/CD4-Analysis/2020-09-22/plots_2020-09-22/cell_type_location_cd4_BCLL-8-T.pdf"
##
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/CD4-Analysis/2020-09-22/plots_2020-09-22/cell_type_location_cd4_BCLL-10-T.pdf"
##
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/CD4-Analysis/2020-09-22/plots_2020-09-22/cell_type_location_cd4_BCLL-12-T.pdf"
##
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/CD4-Analysis/2020-09-22/plots_2020-09-22/cell_type_location_cd4_BCLL-13-T.pdf"
##
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/CD4-Analysis/2020-09-22/plots_2020-09-22/cell_type_location_cd4_BCLL-14-T.pdf"
##
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/CD4-Analysis/2020-09-22/plots_2020-09-22/cell_type_location_cd4_BCLL-9-T.pdf"
##
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/CD4-Analysis/2020-09-22/plots_2020-09-22/cell_type_location_cd4_BCLL-11-T.pdf"
Look at correlation between cell types and markers
sp_sub <- sp_obj[, sp_obj@meta.data$gem_id == sample_id]
sp_sub@images <- sp_sub@images[Seurat::Images(sp_sub) == sample_id]
scatterplot2 <- function(se_obj, x, y, color) {
df <- data.frame(
feat1 = se_obj@assays$MAGIC_Spatial@data[x, ],
feat2 = se_obj@meta.data[, y],
color = se_obj@meta.data[, color]
)
ggplot2::ggplot(df,
ggplot2::aes(x = feat1,
y = feat2)) +
ggplot2::geom_point(ggplot2::aes(color = color)) +
ggplot2::geom_smooth(method = "loess", se = FALSE) +
# ggplot2::geom_smooth(method = lm, se = FALSE) +
ggplot2::theme_classic() +
ggplot2::labs(
x = x,
y = glue::glue("Proportion of\n{y}"),
color = color)
}
bcl6_ls <- lapply(ct, function(ii) {
scatterplot2(
se_obj = sp_sub,
x = "BCL6",
y = ii,
color = "annotation")
})
patchwork::wrap_plots(bcl6_ls, ncol = 5) +
patchwork::plot_layout(guides = "collect") &
ggplot2::theme(legend.position = "right")
prdm1_ls <- lapply(ct, function(ii) {
scatterplot2(
se_obj = sp_sub,
x = "PRDM1",
y = ii,
color = "annotation")
})
patchwork::wrap_plots(prdm1_ls, ncol = 5) +
patchwork::plot_layout(guides = "collect") &
ggplot2::theme(legend.position = "right")
tcf7_ls <- lapply(ct, function(ii) {
scatterplot2(
se_obj = sp_sub,
x = "TCF7",
y = ii,
color = "annotation")
})
patchwork::wrap_plots(tcf7_ls, ncol = 5) +
patchwork::plot_layout(guides = "collect") &
ggplot2::theme(legend.position = "right")
lef1_ls <- lapply(ct, function(ii) {
scatterplot2(
se_obj = sp_sub,
x = "LEF1",
y = ii,
color = "annotation")
})
patchwork::wrap_plots(lef1_ls, ncol = 5) +
patchwork::plot_layout(guides = "collect") &
ggplot2::theme(legend.position = "right")
cxcr5_ls <- lapply(ct, function(ii) {
scatterplot2(
se_obj = sp_sub,
x = "CXCR5",
y = ii,
color = "annotation")
})
patchwork::wrap_plots(cxcr5_ls, ncol = 5) +
patchwork::plot_layout(guides = "collect") &
ggplot2::theme(legend.position = "right")
Seurat::SpatialFeaturePlot(
object = sp_obj,
features = c("CD3D", "CD4", "CD8A", "CD8B"),
alpha = c(0, 1),
images = "esvq52_nluss5",
slot = "data") |
Seurat::SpatialFeaturePlot(
object = sp_obj,
features = c("Cytotoxic", "Naive"),
alpha = c(0, 1),
images = "esvq52_nluss5")